#!/usr/bin/perl
# 2012/12/05 Bug fix. Location was differently reported
# 2012/11/06 Alignment by BWA with several information for HiD or ChIA-PET data
# 2012/10/27 Alignment by BWA without allowing repeat location for HiD (ChIA-PET) data


use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;


#---------------------------------------
# FIX parameter
#---------------------------------------
my $LENGTH = 25;	# length of each reads
#my $FILE_reference = "/home/hideki/bio/bwa/pombe2009";	# reference for BWA
my $FILE_reference = "/home/hideki/bio/bwa/pombe.ASM294v2.19";	# reference for BWA

if(@ARGV != 2 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [sanger fastq file ex.hoge_1.fastq]\n";
}



my %opt;
getopts("i:", \%opt);
my $FILE_fastq1 = $opt{i};
(my $FILE_fastq2 = $FILE_fastq1) =~ s/_1.fastq/_2.fastq/;
my ($name, $dir, $ext) = &fileparse($FILE_fastq1, '\..*');
$name =~ s/_1$//;
my $FILE_sai1 = $dir . $name . '_1.sai';
my $FILE_sai2 = $dir . $name . '_2.sai';
my $FILE_alnLog1 = $dir . $name . '_1_aln.log';
my $FILE_alnLog2 = $dir . $name . '_2_aln.log';
my $FILE_sam1 = $dir . $name . '_1.sam';
my $FILE_sam2 = $dir . $name . '_2.sam';
my $FILE_samLog1 = $dir . $name . '_1_sam.log';
my $FILE_samLog2 = $dir . $name . '_2_sam.log';
my $FILE_map = $dir . $name . '.map';


#---------------------------------------
# aln
#---------------------------------------
unless(-f $FILE_sai1){
	system("bwa aln -t 4 $FILE_reference $FILE_fastq1 > $FILE_sai1 2> $FILE_alnLog1") == 0 or die "aln command failed: $!";
#	system("bwa aln -t 4 -I $FILE_reference $FILE_fastq1 > $FILE_sai1 2> $FILE_alnLog1") == 0 or die "aln command failed: $!";	# only for old illumina
}
unless(-f $FILE_sai2){
	system("bwa aln -t 4 $FILE_reference $FILE_fastq2 > $FILE_sai2 2> $FILE_alnLog2") == 0 or die "aln command failed: $!";
#	system("bwa aln -t 4 -I $FILE_reference $FILE_fastq2 > $FILE_sai2 2> $FILE_alnLog2") == 0 or die "aln command failed: $!";	# only for old illumina
}

#---------------------------------------
# samse
#---------------------------------------
# output repat until 100 locations
unless(-f $FILE_sam1){
	system("bwa samse -n 100 $FILE_reference $FILE_sai1 $FILE_fastq1 > $FILE_sam1 2> $FILE_samLog1") == 0 or die "samse command failed: $!";
}
unless(-f $FILE_sam2){
	system("bwa samse -n 100 $FILE_reference $FILE_sai2 $FILE_fastq2 > $FILE_sam2 2> $FILE_samLog2") == 0 or die "samse command failed: $!";
}


#---------------------------------------
# SAM format
#---------------------------------------
#@SQ	SN:chr1	LN:5579133
#@SQ	SN:chr2	LN:4539804
#@SQ	SN:chr3	LN:2452883
#2_ChIA-PET_linkerAA	0	chr2	1896903	37	25M	*	0	0	ATTAGCCAAGCGGTTTCATGATATC	@DDDBABDFHHEEIIIAHI@9<BHB	XT:A:U	NM:i:0	X0:i:1	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:25
#4_ChIA-PET_linkerAA	16	chr3	2789	0	25M	*	0	0	ATATCCCTCCTCTCTCCACCAAAAG	B??::@FDAFBGHFDFDADA?BB?:	XT:A:R	NM:i:0	X0:i:5	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:25
#5_ChIA-PET_linkerAA	16	chr3	1233976	37	25M	*	0	0	TCCTTTCGAGTATCAACCAGCTTCA	?FCA3AC3IC>BBIFFDF<AA++A1	XT:A:U	NM:i:0	X0:i:1	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:25
#6_ChIA-PET_linkerAA	16	chr3	16411	0	25M	*	0	0	AATATTTTTTTTACTAGGATTTGTG	HBF@IIGIIIIGEIFDFHFDDDDA:	XT:A:R	NM:i:0	X0:i:3	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:25
#7_ChIA-PET_linkerAA	0	chr2	1601563	0	25M	*	0	0	TTTCGATGACTTCCCTAATAAATTA	@FFFDFFFGDHGG>FHIEEHIIII@	XT:A:R	NM:i:0	X0:i:2	X1:i:0	XM:i:0	XO:i:0	XG:i:0	MD:Z:25	XA:Z:chr2,+1646901,25M,0;


#---------------------------------------
# map format
#---------------------------------------
# 1. id
# 2. chr 1
# 3. position 1
# 4. direction 1
# 5. map quality 1
# 6. hindiii id 1
# 7. hindiii location 1
# 8. chr 2
# 9. position 2
# 10. direction 2
# 11. map quality 2
# 12. hindiii id 2
# 13. hindiii location 2
# 14. type U(unique), R(repeat), N(contain N) for left:right
# 15. unique aligned number (1 for unique) left:right
# 16. other alignable locations format "chromosome,+/-location,CIGRmark,missmatch" separated by ";" left:right

#---------------------------------------
# mapping two files
#---------------------------------------
my $TOTAL_read = 0;
my $Aligned_single = 0;
my $Aligned_both = 0;
my $Not_aligned = 0;

my $fh_sam1 = IO::File->new($FILE_sam1) or die "cannot open $FILE_sam1: $!";
my $fh_sam2 = IO::File->new($FILE_sam2) or die "cannot open $FILE_sam2: $!";
my $fh_map = IO::File->new($FILE_map, 'w') or die "cannot write $FILE_map: $!";
while(my $sam1 = $fh_sam1->getline()){
	my $sam2 = $fh_sam2->getline();
	if($sam1 =~ m/^\@/){
		next;
	}
	$sam1 =~ s/\r?\n//;
	$sam2 =~ s/\r?\n//;
	my ($id1, $flag1, $chr1, $position1, $mapQ1, $CIAGR1, $mate_a1, $mate_b1, $mate_c1, $seq1, $quality1, @option1) = split /\t/, $sam1;
	my ($id2, $flag2, $chr2, $position2, $mapQ2, $CIAGR2, $mate_a2, $mate_b2, $mate_c2, $seq2, $quality2, @option2) = split /\t/, $sam2;

	if($id1 ne $id2){
		die ("$id1 and $id2 is different\n");
	}

	$TOTAL_read++;
	my $flag_aligned = 0;
	my ($direction1, $direction2) = ('NA','NA');

	if($flag1 == 0){
		$direction1 = '+';
		$position1 = $position1 + $LENGTH - 1;
		$flag_aligned++;
	}elsif($flag1 == 16){
		$direction1 = '-';
		$flag_aligned++;
	}else{
		$chr1 = 'NA';
		$position1 = 'NA';
	}



	if($flag2 == 0){
		$direction2 = '+';
		$flag_aligned++;
	}elsif($flag2 == 16){
		$direction2 = '-';
		$position2 = $position2 + $LENGTH - 1;
		$flag_aligned++;
	}else{
		$chr2 = 'NA';
		$position2 = 'NA';
	}

	my ($hinID1, $hinLoc1) = ('NA', 'NA');
	my ($hinID2, $hinLoc2) = ('NA', 'NA');

	if($flag_aligned == 0){
		$Not_aligned++;
	}elsif($flag_aligned == 1){
		$Aligned_single++;
	}elsif($flag_aligned == 2){
		$Aligned_both++;
	}

	my ($xt1, $x01, $xa1) = &CheckOption(@option1);
	my ($xt2, $x02, $xa2) = &CheckOption(@option2);

	my $xt = $xt1 . ":" . $xt2;
	my $x0 = $x01 . ":" . $x02;
	my $xa = $xa1 . ":" . $xa2;

	$fh_map->print("$id1\t$chr1\t$position1\t$direction1\t$mapQ1\t$hinID1\t$hinLoc1\t$chr2\t$position2\t$direction2\t$mapQ2\t$hinID2\t$hinLoc2\t$xt\t$x0\t$xa\n");
}
$fh_sam1->close();
$fh_sam2->close();
$fh_map->close();

printf "Single aligned:\t%d\n", $Aligned_single;
printf "Both aligned:\t%d\n", $Aligned_both;
printf "Not aligned:\t%d\n", $Not_aligned;
printf "Total read:\t%d\n", $TOTAL_read;


sub CheckOption{
	my @data = @_;
	my $xt = 'X';
	my $x0 = 0;
	my $xa = 'NA';
	foreach my $d(@data){
		if($d =~ m/^XT:A:(\w)/){
			$xt = $1;
		}elsif($d =~ m/^X0:i:(\d{1,3})/){
			$x0 = $1;
		}elsif($d =~ m/^XA:Z:(.+)/){
			$xa = $1;
		}
	}
	return ($xt, $x0, $xa);
}

